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Abstract 

This paper considers properties of an optimization based sampler for targeting the poste¬ 
rior distribution when the likelihood is intractable. It uses auxiliary statistics to summarize 
information in the data and does not directly evaluate the likelihood associated with the 
specified parametric model. Our reverse sampler approximates the desired posterior distri¬ 
bution by first solving a sequence of simulated minimum distance problems. The solutions 
are then re-weighted by an importance ratio that depends on the prior and the volume of the 
Jacobian matrix. By a change of variable argument, the output are draws from the desired 
posterior distribution. Optimization always results in acceptable draws. Hence when the 
minimum distance problem is not too difficult to solve, combining importance sampling with 
optimization can be much faster than the method of Approximate Bayesian Computation 
that by-passes optimization. 
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THE REVERSE SAMPLER 

1 Introduction 


Maximum likelihood estimation rests on the ability of a researcher to express the joint den¬ 
sity of the data, or the likelihood, as a function of K unknown parameters 6. Inference can 
be conducted using classical distributional theory once the mode of the likelihood function is 
determined by numerical optimization. Bayesian estimation combines the likelihood with a 
prior to form the posterior distribution from which the mean and other quantities of interest 
can be computed. Though the posterior distribution may not always be tractable, it can be 
approximated by Monte Carlo methods provided that the likelihood is available. When the 
likelihood is intractable but there exists L > K auxiliary statistics -0 with model analog 0(0) 
that is analytically tractable, one can still estimate 0 by minimizing the difference between 0 
and 0(0). 

Increasingly, parametric models are so complex that neither the likelihood nor 0(0) is 
tractable. But if the model is easy to simulate, the mapping 0(0) can be approximated by 
simulations. Estimators that exploit this idea can broadly be classihed into two types. One 
is simulated minimum distance estimator (SMD), a frequentist approach that is quite widely 
used in economic analysis. The other is the method of Approximate Bayesian Computation 
that is popular in other disciplines. This method, ABC for short, approximates the posterior 
distribution using auxiliary statistics 0 instead of the full dataset y. It takes draws of 0 from 
a prior distribution and keeps the draws that, when used to simulate the model, produces 
auxiliary statistics that are close to the sample estimates 0. Both the ABC and SMD can be 
regarded as likelihood free estimators in the sense that the likelihood that corresponds to the 
structural model of interest is not directly evaluated. 

While both the SMD and ABC exploit auxiliary statistics to perform likelihood free esti¬ 
mation, there are important differences between them. The SMD solves for the 0 that makes 
0 close to the average of 0(0) over many simulated paths of the data. In contrast, the ABC 
evaluates 0(0) for each draw from the prior and accepts the draw only if 0(0) is close to 0. The 
ABC estimate is the average over the accepted draws, which is the posterior mean. In |Forneron| 


and Ng (2014), we focused on the case of exact identihcation and used a reverse sampler (RS) 


to better understand the difference between the two approaches. The RS approximates the pos¬ 
terior distribution by solving a sequence of SMD problems, each using only one simulated path 


of data. Using stochastic expansions as in Rilstone et al. (1996) and Bao and Ullah (2007), we 
reported that in the special case when 0(0) = 0 (i.e the auxiliary model is the assumed model), 
the SMD has an unambiguous bias advantage over the ABC. But in more general settings, the 
ABC can, by clever choice of prior, eliminate biases that are inherent in the SMD. 

In this paper, we extend the analysis to over-identihed models and provide a deeper un¬ 
derstanding of the reverse sampler. The RS is shown to be an optimization-based importance 
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sampler that transforms the density from draws of to draws of 6 so that when multiplied by 

the prior and properly weighted, the draws follow the desired posterior distribution. Section 2 
considers the exactly identified case and shows that the importance ratio is the determinant of 
the Jacobian matrix. Section 3 considers the over-identified case when the dimension of 
exceeds that of 6. Because of the need to transform densities of different dimensions, the deter¬ 
minant of the Jacobian matrix is replaced by its volume. Using analytically tractable models, 
we show that the RS exactly reproduces the desired posterior distribution. 

The RS was initially developed as a framework to better understand the different approaches 
to likelihood free estimation. While not intended to compete with existing implementations of 
ABC, the use of optimization in RS turns out to have a property that is of independent interest. 
Creating a long sequence of ABC draws such that the simulated statistic -0^ and the data 0 
deviate by no more than 5 can take infinite time if 6 is set to exactly zero as theory suggests. 
This has generated interests within the ABC community to control for 6. The RS by-passes this 
problem because SMD estimation makes 0^ as close to 0^ as machine precision permits. We 
elaborate on this feature in Section 4. Of course, the RS is useful only when the SMD objective 
function is well behaved and easy to optimize, which may not always be the case. But allowing 
optimization to play a role in ABC can be useful, as independent work by Meeds and Welling| 


(2015) also found. 


1.1 Preliminaries 


In what follows, we use a ‘hat’ to denote estimators that correspond to the mode (or extremum 
estimators) and a ‘bar’ for estimators that correspond to the posterior mean. We use (s, S) and 
(b, B) to denote the (specific, total number of) draws in frequentist and Bayesian type analyses 
respectively. A superscript s denotes a specific draw and a subscript S denotes the average 
over S draws. These parameters S and B have different roles. The SMD uses S simulations to 
approximate the mapping 0(0), while the ABC uses B simulations to approximate the posterior 
distribution of the infeasible likelihood. 

We assume that the data y = (yi,..., j/t)^ have finite fourth moments and can be represented 
by a parametric model with probability measure Ve where 0 G 0 C IR^, 0o is the true value. 
The likelihood L(0|y) is intractable. Estimation of 0 is based on L > K auxiliary statistics 
0(y(0o)) which we simply denote by 0 when the context is clear. The model implies statistics 
0(0). The classical minimum distance estimator is 


0CMD = argmin0j(0,0(0)) = g{e)'Wg{e), y(0) = 0 - 0(0). 
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i There exists a unique interior point 6q € Q (compact) that minimizes the popula¬ 
tion objective function — ■j/>(0))'W('j/>(0o) — '0(^))- The mapping 6 —^ ^(0) = 

lim 7 ’_>.oo IE['i/’(0)] is continuously differentiable and injective. The L x K Jacobian matrix 
ipei^) = has full column rank, and the rank is constant in the neighborhood of Oq. 

ii There is an estimator r/) such that y/T{xl} — ■j/>(0o))-^AA(O, S). 

iii W is a L X L positive definite matrix and W'lj^e^Oo) has rank K. 


Assumption A ensures global identification and consistent estimation of 6, see Newey and 


McFadden (1994). In Gourieroux et al. (1993), the mapping ip : 6 ^ '0(^) is referred to as the 


binding function while in Jiang and Turnbull] (2004), ip(6) is referred to as a bridge function. 
When ip{0) is analytically intractable, the simulated minimum distance estimator (SMD) is 


OsMB = argmin0Js(^,^s(0)) = argmin0C?5(0)'Wff5(0). 


( 1 ) 


where S' > 1 is the number of simulations. 


S = 1 


Notably, the term K[ip{0)] in CMD estimation is approximated by ^ '0^(y*(^))- The SMD 

was first used in [Smith] ( |1993 ). Different SMD estimators can be obtained by suitable choice 
of the moments g{0), including the indirect inference estimator of Gourieroux et al. (1993), 


the simulated method of moments of Duffie and Singleton (1993), and the efficient method of 


moments of Gallant and Tauchen (1996). 


The first ABC algorithm was implemented by Tavare et al. (1997) and Pritchard et al 
(1996) to study population genetics. They draw 6^ from the prior distribution simulate 

the model under 0^ to obtain data y^, and accept 0^ if the vector of auxiliary statistics ip{6^) 
deviates from ip by no more than a tuning parameter 6. If ip are sufficient statistics and J = 0, 
the procedure produces samples from the true posterior distribution if B —)> oo. 


The Accept-Reject ABC: For b = 1,..., B 

i Draw from tt{6) and from an assumed distribution Fg, 

ii Generate y^(£^,'i?) and ip^ = ip{y^). 

iii Accept 0^ = r? if — ip'^ W (^ip^ — ip'^ < 6. 
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The accept-reject method (hereafter, AR-ABC) simply keeps those draws from the prior distri¬ 
bution '/r(0) that produce auxiliary statistics which are close to the observed 1 / 7 . As it is not easy 
to choose 6 a priori, it is common in AR-ABC to fix a desired quantile q, repeat the steps [B / 
times. Setting 6 to the g-th quantile of the sequence of that will produce exactly B draws is 


analogous to the idea of keeping fc—nearest neighbors considered in Gao and Hong (2014). 

Since simulating from a non-informative prior distribution is inefficient, the accept-reject 
sampler can be replaced by one that targets at features of the posterior distribution. There 
are many ways to target the posterior distribution. We consider the MCMC implementation of 


ABC proposed in Marjoram et al. (2003) (hereafter, MCMC-ABC). 


The MCMC-ABC: For b = I,..., B with 6^ given and proposal density g( j0^), 

i Generate i? ~ q{'d\0^) 

ii Draw errors from and simulate data i?). Compute -0^^^ = 0(y^^^). 


hi Set to 1 ? with probability PABcif ^^and to 6^^^ with probability 1 — 1 '^) 

where 

( 2 ) 

The AR and MCMC both produce an approximation to the posterior distribution of 0. It is 
common to use the posterior mean of the draws ^ ^ as the ABC estimate. The 

MCMC-ABC uses a proposal distribution to account for features of the data so that it is 
less likely to have proposed values with low posterior probability. The tuning parameter 5 
affects the bias of the estimates. Too small a 5 may require making many draws which can be 
computationally costly. 

The ABC samples from the joint distribution of (0^, 0^(e^, 0^)) and then integrates out 
The posterior distribution is thus 


The indicator function (also the rectangular kernel) equals one if ||0 — 0^|| does not exceed 5. 
The ABC draws are dependent due to the Markov nature of the MCMC-ABC sampler. 

Both the SMD and ABC assume that simulations provide an accurate approximation of 


0(0) and that auxiliary statistics are chosen to permit identification of 0. Creel and Kristensen 


( |2015 ) suggests a cross-validation method for selecting the auxiliary statistics. For the same 
choice of 0, the SMD finds the 0 that makes the average of the simulated auxiliary statistics 
close to 0. The ABC takes the average of 0^, drawn from the prior, with the property that each 
0^ is close to 0. In an attempt to understand this difference, Forneron and Ng (2014), takes 


as starting point that each 0^ in the above ABC algorithm can be reformulated as an SMD 
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problem with S' = 1. We consider an algorithm that solves the SMD problem many times to 

obtain a distribution for 6^, each time using one simulated path. The sampler terminates with 
an evaluation of the prior probability, in contrast to the ABC which starts with a draw from 
the prior distribution. Hence we call our algorithm a reverse sampler (hereafter, RS). The RS 
produces a sequence of 6^ that are independent optimizers and do not have a Markov structure. 

In the next two sections, we explore additional features of the RS. As an overview, the 
distribution of draws that emerge from SMD estimation with S' = 1 may not be from the 
desired posterior distribution. Hence the draws are re-weighted to target the posterior. In the 
exactly identified case, can be made exactly equal to '0 by choosing the SMD estimate as 
6^. Thus the RS is simply an optimization based importance sampler using the determinant of 
Jacobian matrix as importance ratio. In the over-identified case, the volume of the (rectangular) 
Jacobian matrix is used in place of the determinant. Additional weighting is given to those 0° 
that yields '0^ sufficiently close to tp. 


2 The Reverse Sampler: Case K = L 


The algorithm for the case of exact identification is as follows. For b = 1,... ,B 

i Generate from F^. 

ii Find 6^ = argmin^j('('0^(0,e^),■0) and let = 'ip^{0^,e^). 
hi Set u;(6>^£^) = Tr{e^)\'$l{0\ e^)\-\ 
iv Re-weigh the 6^ by 


w{e^) 


Like the ABC, the draws 6^ provides an estimate of the posterior distribution of 6 from which 
an estimate of the posterior mean: 

^ ' Qb\ 


Ors = 


w{0° 


- 0 “ 


^ J2b=iw{o^) 

can be used as an estimate of 0. Each 0^ is a function of the data xp and the draws that 
minimizes j\{xp{0,e^),xp). The K first-order conditions are given by 


F{0\ e\ip) = ^9i{.0\e\ip) ^ ^ 


(3) 


where is the LxK matrix of derivatives with respect to 0 evaluated at the arguments. 

It is assumed that, for all b, this derivative matrix has full column rank K. For SMD estimation, 
~ This Jacobian matrix plays an important role in the RS. 

The importance density denoted h{0^,s^\'ip) is obtained by drawing from the assumed 
distribution and finding 0^ such that .J{xp^{0, e^),'ip) is smaller than a pre-specified tolerance. 
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When K = L, this tolerance can be made arbitrarily small so that up to numerical precision, 

s^) = rp. This density h{0^, e^\ip) is related to p^t, e^)) = p{'tp^, e^) by a change 

of variable: 

Now p{6^^'ip^\'ip) oc p{xp\6^,'ip^)p['ip^^e^\6^)'n{6^) and p{'ip\0^^xp^) is constant since xp^ = xp. 
Hence 

where the weights are, assuming invertibility of the determinant: 


u;(0^e') = 7^(0')|^^(0^£^^)|-^ (4) 

Note that in general, ^ 

In the above, we have used the fact that ]I||^_^i,||_g is 1 with probability one when K = L. 
The Jacobian of the transformation appears in the weights because the draws 6^ are related to 
the likelihood via a change of variable. Hence a crucial aspect of the RS is that it re-weighs the 
draws of 6^ from h{6^,e). Put differently, the unweighted draws will not, in general, follow the 
target posterior distribution. 

Consider a weighted sample (0^, w{9^, s)) with w{6^, e^) defined in Q. The following propo¬ 
sition shows that as H —oo, RS produces the posterior distribution associated with the infea¬ 
sible likelihood, which is also the ABC posterior distribution with J = 0. 


Proposition 1 Suppose thatxp^ : 0 —?■ xp^{9, e^) is one-to-one and the determinant | | , 

\xp^Q{6,e^,xp)\ is hounded away from zero around 6^. For any measurable function p{6) such 
t/iat ((/? (0)) = f (p (0) p(d\xp)d0 exists, then 


Efw(a‘.e‘’)y(#) 


a.s. 


®'p(0|-0) (^)) • 
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Convergence to the target distribution follows from a strong law of large numbers. Fixing 

the event = xj) is crucial to this convergence result. To see why, consider first the numerator: 


JJ (p{e)w{e,e)p{xi^\e’’\e)\xlx0{e,e,'$)\d 




-1 


(p{0) xljf){9,e,xlx) 7r{0)p{xlx ,e \9) xl)(){9,£,xp) 
p (9) tt{9)p{xI}’^, £\9)d£d9 
(p (0) Tr{9)p{xlx, £\9)d£d9 
= I p{9)TT{9)L{xjj\9)d9. 


d£^d9 


Furthermore, the denominator converges to the integrating constant since 
f 7r(9)L(x/xl9)d9. Proposition 0 implies that the weighted average of 9^ converges to the pos¬ 
terior mean. Furthermore, the posterior quantiles produced by the reverse sampler tends to 
those of the infeasible posterior distribution p(0|’j/>) as i? —)• oo. As discussed in 


Forneron and 


Ng (2014), the ABC can be presented as an importance sampler. Hence the accept-reject algo¬ 


rithm in Tavare et al. (1997) and Pritchard et al. (1996), as well as the Sequential Monte-Carlo 


approach to ABC in Sisson et al. (2007); Toni et al. (2009) and Beaumont et al. (2009) are 
all important samplers. The RS differs in that it is optimization based. It is also developed 


independently in Meeds and Welling (2015). 

We now use examples to illustrate how the RS works in the exactly identihed case. 


Example 1: Suppose we have one observation y ~ AA (0,1) or y = 0 -I- e, e ~ W(0,1). The 
prior for 0 is 0 ~ A/" (0,1). By drawing, ~ AA(0,1), we obtain = 6^ + ^ AA(0, 2). The 

ABC keeps 9^\y^ = y. Since {9^,y^) are jointly normal with covariance of 1, we deduce that 
9^\y^ = y M{y/2, 1/2). The exact posterior distribution for 9 is N{y/2, 1/2). 


The RS draws ~ A/" (0,1) and computes 9^ = y — which is M{y, 1) conditional on y. 
The Jacobian of the transformation is 1. Re-weighting according to the prior, we have: 

PRsi0\y) oc 4>i9)4>{9 - y) oc exp (^-1 (6»2 + (6/ _ y)2)^ oc exp (20^ - 26»y)^ 

ocexp (^-‘^{9 -y/2f'^ . 

This is the exact posterior distribution as derived above. 

Example 2 Suppose y = Q{u, 9),e ~ ^o,i] Q is a quantile function that is invertible and 
differentiable in both arguments^ For a single draw, y is a sufficient statistic. The likelihood- 
^We thank Neil Shephard for suggesting the example. 
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based posterior is: 


p{9\y) oc TT{e)f{y\e). 
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The RS simulates y^(9) = Q(e^|0) and sets = y. Or, in terms of the CDF: 

e' = F{y\e^) 

Consider a small perturbation to y holding v!’ fixed: 


Q = dy 


dF{y\e^ 

dy 


,dFiy\e 


+ de^ ' = dyF;^{y\6<’) + 


In the above, / = iT(-) is the density of y given 6. The Jacobian is: 


d9^ 


F'iym 


f{y\0^) 

dy 


F'tivm 


F;,,iym 


To find the distribution of 9^ conditional on y, assume F{y ,.) is increasing in 9: 

= F(Fiy\9^)<F{y\t)\y 


< 


= F\^e° < 

= F{y\t). 

By construction, f{9\y) = Fg(y|0)j^ Putting things togetherj^ 

f{y\o) 


PRs{0\y) oc 7r{9)\Fg{y\9)\ 


= '^{d)fiy\9) (xp{9\y). 


Fl){y\9) 

Example 3: Normal Mean and Variance We now consider an example in which the 


estimators can be derived analytically, and given in Forneron and Ng (2014). We assume 
yt = St ^ N(m,a‘^). The parameters of the model are 6 = We consider the auxiliary 


statistics: V’(y)^ = [y ■ The parameters are exactly identified. 
The MLE of 9 is 


m = - 


t=l 

We consider the prior 7r(m, (j^) = ck > 0 so that the log posterior distribution is 

T 

—T 1 , 

logp(0|m,u^) oc — log(27r)o-^ - alogcr^ - ^ Y^y^ ~ ' 


t=i 


t=i 


^ If Ffy, •) is decreasing in 6, we have < t\y) = 1 — F{y,t). 

^An alternative derivation is to note that t = ¥ {u < t\y) = P (m = F{y, < t\y) ~ I” — F~^{y, t) = t'\y) . 
Hence }{e'’\y) = §r = ^ = F 2 {y,t) as above. 
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Since 'i/’(y) are sufficient statistics, the RS coincides with the likelihood-based Bayesian esti¬ 
mator, denoted B below. This is also the infeasible ABC estimator. We focus discussion on 
estimators for which have more interesting properties. Under a uniform prior, we obtain 


Q _ 

/ ^ _1_ 

b=i i^k=i 

In this example, the RS is also the ABC estimator with 5 = 0. It is straightforward to show that 
the bias reducing prior is a = 1 and coincides with the SMD. Table shows that the estimators 
are asymptotically equivalent but can differ for fixed T. 




'^SMD 


^RS 


Table 1: Properties of the Estimators 

MSE 


Estimator Prior 

m 

Bias 

Variance 

Gml 

9b 1 

9rs 1 

^SMD 

^ T-5 

2r^ 

^ T-5 
,2 S{T-1) 

,72 

T 

20-2 

T-5 

2a2 

T-5 

2 o -2 

r,_4 T—1 

^<2 (' r _ 5)2 

2cj^ "T-l 
(T-5)2 

0^4 1 

^ S{T-l)-2 

5(T-l)-2 

^d\-) 1 17 rjn 


2& 


2T^ 
4 r-n 


2& 


4 T-l-l 
(r-5)2 


A(7^ 

{S{T-1)-2)Y 


where ki( 5, T) = tends to one as S tend to infinity. 


To highlight the role of the Jacobian matrix in the RS, the top panel of Eigure plots 
the exact posterior distribution and the one obtained from the reverse sampler. They are 
indistinguishable. The bottom panel shows an incorrectly constructed reverse sampler that 
does not apply the Jacobian transformation. Notably, the two distributions are not the same. 
Re-weighting by the Jacobian matrix is crucial to targeting the desired posterior distribution. 

Eigure presents the likelihood based posterior distribution, along with the likelihood free 
ones produced by ABC and the RS-JI (just identified) for one draw of the data. The ABC results 
are based on the accept-reject algorithm. The numerical results corroborate with the analytical 
ones: all the posterior distributions are very similar. The RS-JI posterior distribution is very 
close to the exact posterior distribution. Eigure [T] also presents results for the over-identified 
case (denoted RS-OI) using two additional auxiliary statistics; tp = (y, 5^,in 4 /a^) where 
/ifc = lE(y^)- The weight matrix is diag(l, 1,1/2,1/2). The posterior distribution is very close 
to RS-JI obtained for exact identihcation. We now explain how the posterior distribution for 
the over-identified case is obtained. 
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3 The RS: Case L > K: 
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The idea behind the RS is the same when we go from the case of exact to overidentification. The 
precise implementation is as follows. Let (■?/>, be a kernel function and be a tolerance 
level such that ]Ko(^,- 0 '’) = ]I||^_^ 6 ||^o. 

For 6 = 1,..., R 

i Generate from iT. 


ii Find 6^ = argmin^ •0) where -0^ 


hi Set w{6^, e^) = 7 r( 0 ^)vol ( 


b/ab 


-1 


■ 0 )) where: vol('i/> 0 ) = 




e ^0 


iv Re-weigh 6^ by 


w{e'‘) 


Ef=i M0^) ■ 


We now proceed to explain the two changes:- the use of volume in place of determinant in the 
importance ratio, and the need for L — K dimensional kernel smoothing. 

The usual change of variable formula evaluates the absolute value of the determinant of 
the Jacobian matrix when the matrix is square. The determinant then gives the infinitesimal 
dilatation of the volume element in passing from one set of variables to another. The main issue 
in the case of overidentification is that the determinant of a rectangular Jacobian matrix is not 


well defined. However, as shown in Ben-Israel (1999), the determinant can be replaced by the 


volume when transforming from sets of a higher dimension to a lower one|^ For a L x K matrix 
A, its volume, denoted vol( 74 ), is the product of the (non-zero) singular values of A: 


vol(A) = 


j ^J\A'A\ L>K, rank(T) = K 
\^/\AA^\ L<K, rank(74) = L. 


Furthermore, if ^4 = BC, vol(^) = vol(iJ)vol(C'). 

To verify that our target distribution is unaffected by whether we calculate the volume or 
the determinant of the Jacobian matrix when K = L, observe that 




d%l) 


de^' 


(5) 


The K first order conditions defined by ([^ become: 

F{e\ e\ ii,) = e', :0)'W ( ^ 


= 0 . 


( 6 ) 


^From 


Ben-Israel 

(2001 

il 



See also httpTTTwww. encyclopediaofmath.org/index.php/Jacobian 
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Since L = K, W can be set to an identity matrix Ik- Furthermore, = xp since 

Ji{0^) = 0 under exact identification. As ^ is a square matrix when K = L, we can directly 

use the fact that Tg{9^, xp)dO + s^, 'ip)d'ip = 0 to obtain the required determinant: 


• 1^1 = \-Te{e>’,s\xP)-^Tr{e\s\xP)\. 
dip ^ 


(7) 


Now to use the volume result, put A = I^, B = and C = But A is just a A-dimensional 


identity matrix. Hence vol(Iic) = vol 


-1 


vol||^) =vol| 


f-j 

\dxpj' 


dip 

dip 

de 


which evaluates to 


or 


dxp 

-1 

80 

80 


8xp 


which is precisely \xpg{0, as given in ([TjQ Hence in the exactly identified case, there is no 
difference whether one evaluates the determinant or the volume of the Jacobian matrix. 

Next, we turn to the role of the kernel function (■?/>, The joint density h{0^,e^) is 
related to p^,, ^,,{xp{0^, e^)) = p{xp^, e^) through a change a variable now expressed in terms of 
volume: 

h{e, e^\xp) = p{xp\ . vol s\ ^)) 

When L > K, the objective function Wxp — xp^Wi^ = Ji > 0 measures the extent to which xp devi¬ 
ates from xp^ when the objective function at its minimum. Consider the thought experiment that 
= 0 with probability 1, such as enabled by a particular draw of e^. Then the arguments above 
for AT = L would have applied. We would still have = f 7 r( 0 ^)p(’ 0 ^, = 

f 10(0^, £^)h(0^, e^\xp)de^, except that the weights are now defined in terms of volume. Propo¬ 
sition 1 would then extend to the case with L > K. 

But in general 7 ^ 0 almost surely. Nonetheless, we can use only those draws that yield 
Ji{0^) that are sufficiently close to zero. The more draws we make, the tighter this criterion 


can be. Suppose there is a symmetric kernel IK 5 (-) satisfying conditions in Pagan and Ullah 


(1999, p.96) for consistent estimation of conditional moments non-parametrically. Analogous to 
Proposition!^ the volume vol(r/’g(0^,■?/>)) is assumed to be bounded away from zero. Then 
as the number of draws H —)■ 00 , the bandwidth 5{B) —)• 0 and B5{B) —)• 00 with 


ws(B){0^,S^) = t ( 0 *) vo 1 ( ) Ks^B){'ip,^‘’) 


-1 


( 8 ) 


“Using the 


implicit 
Pb 


function theorem to compute the gradient gives the same 
result. Since xp” = ip we have: Tg = —ipg{0’’,e^,ipyWipg{9’’,e‘‘,ip) 

Y.^Ppl,e^{e\e'^)w{Pp-Pp\e\e\Pp)) = -Ppl{e\e\PpywPpl{e\e\Pp). Then = 

vol(Te-i)vol(T^) = vol(V^^(0^e^1A))■Wr^vol(^^(0^e^^))-Vol(^‘(0^e^-0))-l|W| 

vol('0e(0^, Hence the weights are the same when we only consider the draws where Ji = 0 

which are the draws we are interested in. 
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1 

5 


jj ip{e)wo{e,e)vol\^ipl{e,e^;'$)jp{'ijj,e’'\e)dede’’ 

V9(0)7r(6')l||^_^,|I^QVol^^^(0, p(■^^ £'’|0)vol(^^^(0, s^; dOds’’ 

(^(0)7r(0)l||^_^6||=oP('^> ^’'W)dOde^ 

^{e)TT{e)L{'$\e)de. 


Similarly, the integrating constant is consistent as ^ ^ e^)-^ f 7r(0)L('i/^\d)dO. Hence, 

the RS sampler still recovers the posterior distribution with the infeasible likelihood. Note that 
the kernel function was introduced for developing a result analogous to Proposition 1, but no 
kernel smoothing is required in practical implementation. What is needed for the RS in the over¬ 
identified case is B draws with sufficiently small Hence, we can borrow the idea used in 

the AR-ABC. Specihcally, we fix a quantile q, repeat [B/q] times until the desired number of 
draws is obtained. Discarding some draws seems necessary in many ABC implementations. 

In summary, there are two changes in implementation of the RS in the over-identified case: 
the volume and the kernel function. Kernel smoothing has no role in the RS when K = L. It is 
interesting to note that while the ABC and RS both rely on the kernel to keep draws close 
to "0^ in the over-identified case, the non-parametric rate at which the sum converges to the 
integral are different. The RS uses the first order conditions e^)'W e^) — = 0 

to indicate which K combinations of -0^(0^, — xp are set to zero, rendering the dimension of 

the smoothing problem L — K. To see this, note first that each draw 0^ from the RS is consistent 


for 00 and asymptotically normal as shown in Forneron and Ng (2014). In consequence, the first 

'ilVw (xp\0\e^)-xp^ = 0, 


/ dip{e) I ^ /_ 

d0 \e=eo~'~^P^VT^J 


order condition (FOC) can be re-written as: ^ 


or 


dO '^=^0 




Since l6>=go^ rank, there exists a subspace of dimension K such that xp^{0^,e^) —xp 

is zero asymptotically. Hence the kernel smoothing problem is effectively L — K dimensional. 
The ABC does not use the FOC. Even in the exactly identified case, the kernel smoothing is 
a L = K dimensional problem. In general, the convergence rate of the ABC is L > AT, the 
dimension of xp. 

The following two examples illustrate the properties of the ABC and RS posterior distri¬ 
butions. The first example uses sufficient statistics and the second example does not. Both 
the ABC and RS achieve the desired number of draws by setting the quantile, as discussed in 
Section 2. 
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Example 4: Exponential Distribution Let yi,... ,yT ^ £{0),T = 5, 0o = 1/2- Now ^ = y 

is a sufficient statistic for yi,..., yx- For a flat prior tt{6) oc le>o we have: 

p{6\y) ocp((9|7/i,.. .,yT) = exp(-6'^y) ~ r(r + l,Ty) 


In the just identified case, we let ^o,l] and yt=— log(l — u\)/9^. This gives: 

r rji / J Wt rp / ^ 


t=l 


t=l 


eb 


Since y^ = y, the Jacobian matrix is: 


Moi = 




dd 


T 

E 


1 log(l - u\) 


^ t=i 




Hence for a given T, the weights are: w{6^,u^) oc = y. We verified that the numerical 

results agree with this analytical result. 

In the over identified case, we consider two moments: 


-0 


b 




b\2 


ELM 




T 

t=l 



Since ^ 


log(l--»t) 

( e '>)2 


— 1^. If <5 = 0, the Jacobian matrix is 


^9 


( 


IspT yl 
T 9 


El 

ye*’ T 


b\2 


ELM) 


2^ 



The volume to be computed is vol('0g) = y li/jg'-j/jgl, as stated in the algorithm. Even if W = I, 
the volume is the determinant of in the exactly identified case, plus a term relating to the 
variance of y^. We computed 'il^g for draws with Jj* ~ 0 using numerical differentiatioij^ and 
verified that the values are very close to the ones computed analytically for this example. 

Figure depicts a particular draw of the ABC posterior distribution (which coincides with 
the likelihood-based posterior since the statistics are sufficient), along with two generated by 
the RS sampler. The first one uses the sample mean as auxiliary statistic and hence is exactly 
identified. The second uses two auxiliary statistics: the sample mean and the sample variance. 
For the AR-ABC, we draw from the prior ten million times and keep the ten thousand nearest 
draws. This corresponds to a value of J = 0.0135. For the RS, we draw one million time^and 

®In practice, since the mapping 6 —>■ ■0^(0) is not known analytically, the derivatives are approximated using 
finite differences: dg^ip’’{6) ~ ^ {s+e^e)^^ (g e^e) g ~ Q. 

^This means that we solve the optimization problem one million times. Given that the optimization problem 
is one dimensional, the one dimensional R optimization routine optimize is used. It performs a combination of 
the golden section with parabolic interpolations. The optimum is found, up to a given tolerance level (the default 
is 10“"^), over the interval [0,10]. 
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keep the ten thousand nearest draws which corresponds to a (5 = 0.0001. As for the weight matrix 

W, if we put Wii > 0 and zero elsewhere, we will recover the exactly identified distribution. 
Here, we intentionally put a positive weight on the variance (which is not a sufficient statistic) 
to check the effect on the posterior mean. With Wn = 1/5 and W 22 = 4/5, the RS posterior 
means are 0.7452 and 0.7456 for the just and overidentifed cases. The corresponding values are 
are 0.7456 and .7474 for the exact posterior and the ABC-AR. They are very similar. 


Example 5: ARMA(1,1): For t = 1, ..., T = 200 and Oq = (ao, Oq, ctq) = (0.5,0.5,1.0), the 
data are generated as 

yt = ayt-i + £t + Ost-i, ej ~ AA(0, a^). 

Least squares estimation of the auxiliary model 


yt — + <i>2yt-2 + (t'syt-s + (k^yt-A + ut 


yields L = 5 > K = 3 auxiliary parameters 


"0 — {4>1,4>2,4>3,4>A,^1)- 


We let 7r{a,6,a) = ]lQ,_ 0 g[_i^i]^o->o W = I 5 which is inefficient. In this example, 0 are not 
sufficient statistics since yt has an infinite order autoregressive representation. 

We draw a from a uniform distribution on [0,3] since ^o,oo] is not a proper density. The 
weights of the RS are obtained by numerical differentiation. The likelihood based posterior 
is computed by MCMC using the Kalman Filter with initial condition eg = 0. As mentioned 
above, the desired number of draws is obtained by setting the quantile instead of setting the 
tolerance 6. For the RS, we keep the 1/10=10% closest draws corresponding to a <5 = 0.0007. 
The Sequential Monte-Carlo implementation of ABC (SMC-ABC) is more efficient at targeting 
the posterior than the ABC-AR. Hence we also compare the RS with SMC-ABC as implemented 
in the Easy-ABC package of Lenormand et al. (2013) ^ The requirement for 10,000 posterior 
draws are as follows: 



AR-ABC 

SMC-ABC 

RS 

Likelihood 

Computation Time (hours) 

63 

25 

5 

0.1 

Effective number of draws 

100 , 000,000 

36,805,000 

10,153,108 


6 

0.0132 

0.0283 

0.0007 



°We implemented the SMC-ABC in two ways. First, we use the procedure irVo et al. (20151 using code 
generously provided by Christopher Drovandi. We also use the Easy-ABC package in R of [Lenormand et al.| 
(20131. We thank an anonymous referee for this suggestion. 
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The difference, both in terms of computation time and number of model simulations, is notable. 

As shown in figure the quality of the approximation is also different, especially for a and a. 
The difference can be traced to b. The b used for the SMC-ABC is effectively much larger than 
for the RS. A better approximation requires a smaller b which implies longer computational 
time. Alternatively stated, the acceptance rate at a low value of b is very low. The caveat is 
that the speed gain is possible only if the optimization problem can be solved in a few iterations 
and reasonably fast. In practice, there will be a trade-off between the number of draws and the 
number of iterations in the optimization step as we further explore below. 


4 Acceptance Rate 


The RS was initially developed in Forneron and Ng (2014) as a framework to help understand 
frequentist (SMD) and the Bayesian (ABC) way of likelihood-free estimation. But it turns out 
that the RS has one computation advantage that is worth highlighting. The issue pertains to 
the low acceptance rate of the ABC. 

As noted above, the ABC exactly recovers the posterior distribution associated with the 
infeasible likelihood if -0 are sufficient statistics and <5 = 0 as noted in 


Blum 


( |2010[ ). Of course, 

(5 = 0 is an event of measure zero, and the ABC has an approximation bias that depends on b. In 
theory, a small b is desired. The ABC needs a large number of draws to accurately approximate 
the posterior and can be computationally costly. 

To illustrate this point, consider estimating the mean m in Example 3 with = \ assumed 
to be known, and 7r(m) oc 1. All computations are based on the software package R. From a 
previous draw a random walk step gives m* = -|- e, e ~ A/'(0,1). For small b, we can 
assume m*\m ~ N{fh, l/T). From a simulated sample of T observations, we get an estimated 
mean fh* ~ l/T). As is typical of MCMC chains, these draws are serially correlated. To 

see that the algorithm can be stuck for a long time if m* is far from m, observe that the event 
m* £ [fh — b, fh + b] occurs with probability 


m* E 


[fh—b, m-|-(5]) = $ (^VT{fh -|- <5 — m*)^ —<I> (^VT{m — b — 



The acceptance probability ¥{fh* E [m — <5, m -|- b])dm* is thus approximately linear in b. 
To keep the number of accepted draws constant, we need to increase the number of draws as 
we decrease b. 

This result that the acceptance rate is linear in b also applies in the general case. Assume 
that xl}*{0*) ~ S/T). We keep the draw if ||0 — 0*(0*)|| < b. The probability of this 

event can be bounded above by Ylf=i i-e- 
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The acceptance probability is still at best linear in 6. In general we need to increase the number 

of draws at least as much as 5 declines to keep the number of accepted draws fixed. 


Table 2: Acceptance Probability as a function of 6 


5 

10 

1 

0.1 

0.01 

0.001 


0.72171 

0.16876 

0.00182 

0.00002 

<0.00001 


Table shows the acceptance rate for Example 3 for Oq = ~ (0,2), T = 20, and 

weighting matrix W = diag(a^, 2u^)/T, 7r(m, cr^) oc ]lo-2>o- results confirm that for small 
values of 6, the acceptance rate is approximately linear in 6. Even though in theory, the targeted 
ABC posterior should be closer to the true posterior when 6 is small, this may not be true in 
practice because of the poor properties of the MCMC chain. At least for this example, the 
MCMC chain with moderate value of 6 provides a better approximation to the true posterior 
density. 


To overcome the low acceptance rate issue, Beaumont et al. (2002) suggests to use local re¬ 


gression techniques to approximate <5 = 0 without setting it equal to zero. The convergence rate 


is then non-parametric. Gao and Hong (2014) analyzes the estimator of Creel and Kristensen 


(2013) and finds that to compensate for the large variance associated with the kernel smoothing. 


the number of simulations need to be larger than to achieve \/T convergence, where K 

is the number of regressors. Other methods that aim to increase the acceptance rate include 


the ABC-SMC algorithm of Sisson et al. (2007); Sisson and Fan (2011), as well as the adaptive 


weighting variant due to Bonassi and West (2015), referred to below as SMC-AW. These meth¬ 


ods build a sequence of proposals to more efficiently target the posterior. The acceptance rate 
still declines rapidly with <5, however. 

The RS circumvents this problem because each 6^ is accepted by virtue of being the solution 
of an optimization problem, and hence — is the smallest possible. In fact, in the exactly 

identified case, 5 = = 0. Furthermore, the sequence of optimizers are independent, and the 

sampler cannot be stuck. We use two more examples to highlight this feature. 


Example 6: Mixture Distribution Consider the example in Sisson et al. (2007), also 


considered in Bonassi and West (2015). Let 'k{9) oc l0g[_io,io] 


x\e ~ l/2A7(6', 1) -L 1/100) 

Suppose we observe one draw x = 0. Then the true posterior is 9\x ~ 1/2A7(0,1)-|-1/2A7(0,1/100) 
truncated to [—10,10]. As in Sisson et al. (2007) and [Bonassi and West (2015), we choose three 
tolerance levels: (2,0.5, 0.025) for AR-ABC. Figure shows that the ABC posterior distribu¬ 
tions computed using accept-reject sampling with 6 = 0.025 are similar to the ones using SMC 
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with and without adaptive weighting. The RS posterior distribution is close to both ABC-SMC 


and ABC-SMC-AW, and all similar to Figure 3 reported in Bonassi and West (2015). However, 
they are quite different from the AR-ABC with 6 = 2 and 0.5 are 2 , showing that the choice of 
6 is important in ABC. 

While the SMC, RS, and ABC-AR sampling schemes can produce similar posterior dis¬ 
tributions, Table shows that their computational time differ dramatically. The two SMC 
algorithms need to sample from a multinomial distribution which are evidently more time con¬ 
suming. When 6 = 0.25, the AR-ABC posterior distribution is close to the ones produced 
by the SMC samplers and the RS, but the computational cost is still high. The AR-ABC is 
computationally efficient when 6 is large, but as seen from Figure the posterior distribution 
is quite poorly approximated. The RS takes 0.0017 seconds to solve, which is amazingly fast 
because for this example, the solution is available analytically. No optimization is involved, and 
there is no need to evaluate the Jacobian because the model is linear. Of course, in cases when 
the SMD problem is numerically challenging, numerical optimization can be time consuming 
as well. Our results nonetheless suggest a role for optimization in Bayesian computation; they 
need not be mutually exclusive. Combining the ideas is an interesting topic for future research. 


Table 3: Computation Time (in seconds) 


RS 

ABC-AR 

ABC-SMC 


J =2 

<5=.5 

J=.025 

Sisson et-al 

Bonassi-West 

.0017 

0.4973 

1.6353 

33.8136 

190.1510 

199.1510 


Example 7: Precautionary Savings The foregoing examples are simple and are serve 
illustrative purposes. We now consider an example that indeed has an infeasible likelihood. In 


Deaton (1991), agents maximize expected utility Eq subject to the constraint 


that assets at+i = (1 -|- r){at + yt — ct) are bounded below by zero, where r is interest rate, 
y is income and c consumption. The desire for precautionary saving interacts with borrowing 
constraints to generate a policy function that is not everywhere concave, but is a piecewise linear 
when cash-on-hand is below an endogenous threshold. The policy function can only be solved 
numerically at assumed parameter values. SMD estimation thus consists of solving the model 


and simulating S auxiliary statistics at each guess 6. Michaelides and Ng (2000) evaluate the 
finite sample properties of several SMD estimators using a model with similar features. Since the 
likelihood for this model is not available analytically. Hence Bayesian estimation of this model 
has not been implemented. Here, we use the RS to approximate the posterior distribution. 

yt ~ iid AA(|U, (T^) with 




We generate T = 400 observations assuming that U{c) = 
r = 0.05, /3 = 10/11, /i = 100, cj = 10 ,7 = 2 as true values. We estimate 6 = ( 7 , /r,(T) and 

















19 


THE REVERSE SAMPLER 

assume (/3, r) are known. We use 10 auxiliary statistics: 

^={y fj,y(0) f,,(0) ree(O) f,,(l) f„,(l) fec(2) f,,(2) f,y{0) f,y(0))' 

where ra;,(j) = ^ ~ — &)• We generate B = 13,423 draws and keep the 3, 356 

(25%) nearest draws to rp. After weighting using the volume of the Jacobian matrix we have 
an effective sample size of 1,421 drawsj^ We use an identity weighting matrix so = 

g{0y~g[6). The Jacobian is computed using finite differences for the RS. As benchmark, we also 
compute an SMD with S = 100, Js = gsi^Ydsi^)- this exercise, the SMD only needs to 
solve for the policy function once at each step of the optimization. Hence the binding function 
can be approximated using simulated data at a low cost. For this example, the programs are 
coded in python. The Nelder-Mead method is used for optimization. 

Table 4: Deaton Model: RS, SMD with W = I 



Posterior Mean/Estimate 

Posterior SD/SE 


7 

a 

7 


a 

RS 

1.86 99.92 

10.48 

0.19 

0.84 

0.37 

SMD 

1.76 99.38 

10.31 

0.12 

0.60 

0.34 


Figure shows the posterior distribution of the RS (blue) along with the SMD distribution 
(purple) as approximated by W(0 smD) Vsmd/^) according to asymptotic theory. Tablej^shows 
that the two sets of point estimates are similar. As explained in Forneron and Ng ( |2014[ ), the 
SMD uses simulations to approximate the binding function while the RS (and by implication 
the ABC) uses simulations to approximate the infeasible posterior distribution. In this example, 
the difference in bias is quite small. We should note that the RS took well over a day to solve 
while the SMD took less than three hours to compute. Whether we use our own code for the 
ABC-MCMC or from packages available, the acceptance rate is too low for the exercise to be 
feasible. 


5 Conclusion 


This paper studies properties of the reverse sampler considered in Forneron and Ng (2014) for 


likelihood-free estimation. The sampler produce draws from the infeasible posterior distribution 
by solving a sequence of frequentist SMD problems. We showed that the reverse sampler uses 
the Jacobian matrix as importance ratio. In the over-identified case, the importance ratio can 
be computed using the volume of the Jacobian matrix. The reverse sampler does not suffer 
from the problem of low acceptance rate that makes the ABC computationally demanding. 


®The effective sample size is computed as '^b where the weights satisfy w’*’ = 1. 
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Figure 1; Normally Distributed data 
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Figure 2: The Importance Weights in RS 
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Figure 3: Exponential Distribution 
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Figure 4: ARMA Model 
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Figure 5; Mixture Distribution 
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Figure 6: Deaton Model: RS and SMD 
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Note: Blue density; RS posterior, Black line: large sample approximation for the SMD estimator 
(identity weighting matrix). 















